% Monomials_2.m is a routine that constructs integration nodes and weights 
% under N-dimensional monomial (non-product) integration rule with 2N^2+1 
% nodes; see supplement to "Numerically Stable and Accurate Stochastic 
% Simulation Approaches for Solving Dynamic Economic Models" by Kenneth L. 
% Judd, Lilia Maliar and Serguei Maliar, (2011), Quantitative Economics 2/2, 
% 173–210 (henceforth, JMM, 2011).
%
% This version: July 14, 2011. First version: August 27, 2009. 
% -------------------------------------------------------------------------
% Inputs:  "N" is the number of random variables; N>=1;
%          "vcv" is the variance-covariance matrix; N-by-N;

% Outputs: "n_nodes" is the total number of integration nodes; 2*N^2+1;
%          "epsi_nodes" are the integration nodes; n_nodes-by-N;
%          "weight_nodes" are the integration weights; n_nodes-by-1
% -------------------------------------------------------------------------
% Copyright © 2011 by Lilia Maliar and Serguei Maliar. All rights reserved. 
% The code may be used, modified and redistributed under the terms provided 
% in the file "License_Agreement.txt".
% -------------------------------------------------------------------------


function [n_nodes,epsi_nodes,weight_nodes] = Monomials_2(N,vcv)

n_nodes = 2*N^2+1;    % Total number of integration nodes

% 1. N-dimensional integration nodes for N uncorrelated random variables with 
% zero mean and unit variance
% ---------------------------------------------------------------------------   

% 1.1 The origin point
% --------------------
z0 = zeros(1,N);       % A supplementary matrix for integration nodes: the 
                       % origin point 

% 1.2 Deviations in one dimension
% -------------------------------
z1 = zeros(2*N,N);     % A supplementary matrix for integration nodes; 
                       % n_nodes-by-N
                       
for i = 1:N            % In each node, random variable i takes value either
                       % 1 or -1, and all other variables take value 0
    z1(2*(i-1)+1:2*i,i) = [1; -1];  
end                    % For example, for N = 2, z1 = [1 0; -1 0; 0 1; 0 -1]

% 1.3 Deviations in two dimensions
% --------------------------------
z2 = zeros(2*N*(N-1),N);   % A supplementary matrix for integration nodes; 
                           % 2N(N-1)-by-N

i=0;                       % In each node, a pair of random variables (p,q)
                           % takes either values (1,1) or (1,-1) or (-1,1) or    
                           % (-1,-1), and all other variables take value 0 
for p = 1:N-1                           
    for q = p+1:N      
        i=i+1;               
        z2(4*(i-1)+1:4*i,p) = [1;-1;1;-1];    
        z2(4*(i-1)+1:4*i,q) = [1;1;-1;-1];
    end
end                 % For example, for N = 2, z2 = [1 1;1 -1;-1 1;-1 1]

% z = [z0;z1*sqrt(N+2);z2*sqrt((N+2)/2)];   % Integration nodes 
                       
% 2. N-dimensional integration nodes and weights for N correlated random 
% variables with zero mean and variance-covaraince matrix vcv 
% ----------------------------------------------------------------------  
sqrt_vcv = chol(vcv);            % Cholesky decomposition of the variance-
                                 % covariance matrix
                                 
R = sqrt(N+2)*sqrt_vcv;          % Variable R; see condition (B.8) in the  
                                 % Supplement to JMM (2011)
                                 
S = sqrt((N+2)/2)* sqrt_vcv;     % Variable S; see condition (B.8) in the  
                                 % Supplement to JMM (2011)
                                 
epsi_nodes = [z0;z1*R;z2*S]; 
                                 % Integration nodes; see condition (B.8)   
                                 % in the Supplement to JMM (2011); 
                                 % n_nodes-by-N

% 3. Integration weights
%-----------------------
weight_nodes = [2/(N+2)*ones(size(z0,1),1);(4-N)/2/(N+2)^2*ones(size(z1,1),1);1/(N+2)^2*ones(size(z2,1),1)];
                                 % See condition in (B.8) in the Supplement  
                                 % to JMM (2011); n_nodes-by-1; the weights 
                                 % are the same for the cases of correlated  
                                 % and uncorrelated random variables